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We have derived equations for nonadiabatic Ehrenfest molecular dynamics which conserve the total energy 
in the case of time-dependent discretization for electrons. A discretization is time-dependent in all cases 
where it or part of it depends on the positions of the nuclei, for example, in atomic orbital basis sets, and 
in the projector augmented-wave (PAW) method, where the augmentation functions depend on the nuclear 
positions. We have derived, implemented, and analyzed the energy conserving equations and their most 
common approximations for a ID test system where we can achieve numerical results converged to a high 
accuracy. Based on the observations in ID, we implement and analyze the Ehrenfest molecular dynamics 
in 3D using the PAW method and the time-dependent density functional formalism. We demonstrate the 
applicability of our method by carrying out calculations for small and medium sized molecules in both the 
adiabatic and the nonadiabatic regime. 
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I. INTRODUCTION 

Many natural processes, such as light absorption, igni- 
tion of chemical reactions and ion-atom collisions, are re- 
lated to excited electronic states and their time develop- 
ment. The most general approach for treating such nona- 
diabatic processes, which in general involve two or more 
coupled electronic states, would be to solve the time- 
dependent many-body Schrodinger equation. However, 
this is not feasible for systems consisting of more than a 
few electrons. Moreover, the standard ab initio molecular 
dynamics (AIMD) methods such as Car-Parrinello MD^ 
(CPMD) or Born-Oppenheimer MD (BOMD) (see, for 
example, Ref. l2) confine electrons to a single adiabatic 
state, typically the ground state. Semiclassical meth- 
ods, such as Ehrenfest molecular dynamics (Ehrenfest 
MD) or trajectory surface hopping (TSH), in which the 
electrons are treated quantum-mechanically via the time- 
dependent Schrodinger equation and the nuclei classi- 
cally via the Newtonian mechanics, have been developed 
some decades ago, but only during the last decade they 
have become feasible in atomistic simulations beyond 
small systems. This is partly due to the methodological 
advances in time-dependent density functional theory'^ 
(TDDFT) which provides a computationally affordable 
basis for the Ehrenfest MD and the TSH methods, and 
also due to the rapidly increasing amount of computa- 
tional resources available. 

Ehrenfest MD within TDDFT offers a simple yet ef- 
fective framework for simulating nonadiabatic processes 
by coupling the time-dependent Kohn-Sham (TDKS) 
equations'* with classical equations of motion for nuclei 
via the KS potential energy surface (PES). The method 
works well for condensed matter, where many single elec- 
tron levels are involved and a single reaction path dom- 
inates the nonadiabatic process such as in carbon nano- 
structures. However, when reactions pass regions of close 
lying electronic states but end up in a state which is well 



described by a single potential energy surface, the TSH 
method is preferable^. This is due to the deficiency of the 
Ehrenfest MD that the system remains in a mixed state 
after exiting the nonadiabatic region. Moreover, due to 
its mean-field character, the Ehrenfest MD cannot cor- 
rectly describe multiple reaction paths^. 

The Ehrenfest MD has been succesfully used for study- 
ing various nonadiabatic processes such as collisions bet- 
ween atomic oxygen and graphite clusters^, excited car- 
rier dynamics in carbon nanotubesS and electronic exci- 
tations in ion bombardment of carbon nanostructures^. 
Previously, there have been only a few implementations 
of the Ehrenfest MD, all of which are based on pseu- 
dopotentials. Various basis sets such as LCAOlSj plane 
wavesii and real space grids^^ have been used. Real space 
techniques have several advantages: 1) the quality of the 
results can be controlled by a single variable, the grid 
spacing, 2) parallelization can be done straightforwardly 
using domain decomposition and 3) different boundary 
conditions such as Dirichlet, periodic, or a mixture of 
them can be easily applied. 

In spite of the extensive use of the the projector 
augmented-wavei^ method in electronic structure calcu- 
lations, to our knowledge, it has not been used in Ehren- 
fest MD simulations previously. Compared to pseudopo- 
tentials, the PAW method improves the description of 
the transition metal elements and the first row elements 
with open p-shells. Moreover, in order to perform ac- 
curate calculations, the PAW method in real space al- 
lows one to use fewer grid points than pseudopotentialsi^ 
as well as longer time steps in the propagation of the 
TDKS equationSii^. The all-electron (AE) nature of 
the PAW method, albeit with the core states frozen, is 
also a methodological advantage over pseudopotentials. 
However, compared to pseudopotentials, Ehrenfest MD 
within the PAW method is more complicated due to the 
augmentation functions that depend on the atomic po- 
sitions. First, an additional term describing the moving 
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spatial gauge of the electrons must be included in the 
TDKS equations—. Also, the Hellmann-Feynman (HF) 
theorem (see, for example, Ref. [13), traditionally used 
for calculating the atomic forces, is no longer valid. 

The present paper is constructed as follows. In Sec. 
we investigate various expressions for the Ehrcnfcst MD 
forces and derive results regarding the total energy con- 
servation in the quantum-classical dynamics. Moreover, 
we describe our propagation algorithm for the Ehren- 
fest MD equations. In Sec. IIII[ we present the one- 
dimensional test system and the example calculations 
carried out with it. The Ehrenfest MD implementation 
within the PAW method is described in Sec. IIVI The 
applicability of our method is demonstrated by carry- 
ing out simulations for the NaCl dimer and the II2C = 
NH^ and C40H16 molecules. Furthermore, the applica- 
bility of the different forces in both the adiabatic and the 
nonadiabatic regime is discussed. Finally, we give brief 
conclusive remarks in Sec. |Vl 



II. METHODOLOGY 

Ehrenfest MD in general can be defined by the time- 
dependent Schrodinger equation for the electrons and the 
Newton equations of motion for the nuclei (atomic units 
are used throughout this paper) 

M,R, = -Vr^E,s = -Vr„ (^-iHeil*) , (2) 

where 5* is the many-particle electronic wavefunction, de- 
pending explicitly on the time t and the electronic degrees 
of freedom {rj}, and implicitly on the atomic positions 
{Ra}. Hei is the electronic Hamiltonian. Thus, the force 
on the nuclei is calculated as an average over all elec- 
tronic adiabatic states, i.e., the Ehrenfest MD is effec- 
tively a mean-field theory. By expanding the electronic 
wavefunctions in the basis obtained by solving the time- 
independent Schrodinger equation, one can show that 
nonadiabatic effects are included in the Ehrenfest MD^. 
In this section, we investigate the Ehrenfest MD within 
the single-particle formalism in a finite basis. 



A. A general time-dependent quantum-classical system 



where is the kinetic energy of the non-interacting elec- 
trons, Eext is the energy due to the external potential 
which we assume to depend explicitly on the atomic po- 
sitions, and E' is a functional that contains nonlinear 
density-dependent terms. The electronic density p is a 
function of atomic positions {Ra} and time t. We use 
the semicolon in Eq. @ to distinguish between the 
function and vector dependencies of E^i- In the Kohn- 
Sham formalism, E' contains the Hartree and exchange- 
correlation energies. H' [Eq. ([3])] is then the functional 
derivative of E'[p] with respect to the density. We ex- 
pand the single-particle electronic states on a basis {xk } 
that depends explicitly on the atomic positions. 



V'„(r, {Ra},t) = Vc„fe(t)xfe(r, {Ra}). 



k 



(5) 



Moreover, we assume that the time-dependency of the ba- 
sis functions is solely due to the movement of the atomic 
positions. With this construction, we define the Hamil- 
tonian matrices H. H° and H' as follows 



H[,^{X^H'\X,)- 



(6) 
(7) 
(8) 



We write the total Hamiltonian of the quantum- 
classical system consisting of quantum-mechanical elec- 
trons and classical nuclei as 



n = 



Pa 



^ 2M„, 



E,: 



(9) 



In quantum-classical molecular dynamics, it is essential 
that the time derivative of the total Hamiltonian % is at 
least approximatively zero, i.e., H, is an invariant. The 
time derivative of the total Hamiltonian reads as 



dU 



dE^ 
dt ' 



(10) 



where Va is the velocity of atom a, and Fq is the atomic 
force. Now, the time derivative of the electronic energy 
[Eq. Q] can be computed as 



dE,i 
dt 



dE,i dE, 



dRa 



dt 



(11) 



We consider a general time-dependent electronic sys- 
tem within Kohn Sham-like single-particle formalism 
when the Hamiltonian operator H can be written as a 
sum of a position-dependent term H° and a nonlinear 
term H' that contains terms depending on the electronic 
density p 



i?(r, {Ra}, p) = i?"(r, {Ra}) + i?'(r, p). 



(3) 



The energy functional of the system is defined in terms 
of p as 

E,i[p; {Ra}] = T,[p] + E,^t[p; {Ra}] + E'[p], (4) 



In order to proceed from the above expression, we expand 
the electronic states tpn in the time-dependent single- 
particle Schrodinger equation. 



(12) 



on the basis {xk} [Eq- ©]. Consequently, we arrive at 
the following matrix representation of the single-particle 
TD Schrodinger equation [Eq. (IT^ ] 



iS 



dt 



(H + P)c„ 



(13) 
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where the vector c„ contains the basis function coefH- 
cients of state n, c„fc, and the matrices S and P are 
defined as 



(14) 
(15) 



The S matrix describes the overlap between the basis 
functions, while the P matrix takes into account the mov- 
ing spatial gauge of the electrons due to the changing 
atomic positions and conserves the norm of the electronic 
states. It disappears if the nuclei do not move. Further- 
more, employing the chain rule, the P matrix can be 
written in terms of the atomic velocities as 



D„ 



with the definition 



) 



(16) 



(17) 



Using the matrix representation of the TD Schrodinger 
equation [Eq. (fT3|)]. we can compute the time derivatives 
of the vectors c„. Moreover, as the nonlinear energy term 
E' [Eq. (lU] is essentially of the form 



E' = j F{v,p{v))dv, 



(18) 



where F is a, density-dependent function. Using the chain 
rule, we can compute the (partial) time derivative of the 
nonlinear energy. 



dE' 
~dt 



(19) 



Since dF/ dp equals the functional derivative of E' with 
respect to the density, Eq. can be written in the 

following matrix form 



dE' 
~dt 



dt 



C.C.]. 



(20) 



Similar reasoning can be applied to the linear part of 
the electronic energy, = Ec\ — E' . Thus, we get the 
following expression for its time derivative 



dE° 



- C.C.]. 



(21) 



With Eqs. ([20| and (1211) and the matrix representation of 
the single-particle TD Schrodinger equation [Eq. 
the time derivative of the electronic energy [Eq. 
takes the following form 



m] 



dE^ 
dt 



dE,i dE, 
dRa dt 



dE,i 
dRa 



E 



<(HS-iD,+c.c.)c„ 



(22) 



The dynamics of the quantum-classical system is defined 
by the Newton equations for the nuclei and Eq. (|13p for 
the electrons, i.e., 



MaRa 



iS 



dCn 

dt 



--Fa, 

= (H 



P)c„. 



(23) 
(24) 



Now, whether the total Hamiltonian usually inter- 
preted as the total energy of the system, is a conserved 
quantity or not, depends on the type of force used in 
the calculations. Next, we will consider three different 
forces: 1) Total energy conserving (EC) force, 2) Incom- 
plete basis set corrected (IBSC) force and 3) Hellmann- 
Feynman (HE) force. 



1. The total energy conserving force 

The most general Ehrenfest MD force expression can 
be derived using the requirement that the time derivative 
of the total Hamiltonian [Eq. ([TUl) ] equals zero. This can 
be achieved by defining the force by the terms inside the 
brackets in Eq. ([22|) with an opposite sign, i.e.. 



dEA 
dRa 



<(HS-iDa + c.c.)c„ 



(25) 



A very similar expression has been derived using first 
the HE theorem to find an equation for the atomic forces 
according to the Ehrenfest MD, and then writing the 
resulting equations in a basis^-. Moreover, Di Ventrai^ 
suggested that for time-dependent processes, one should 
calculate the force acting on atom a as the time derivative 
of the expectation value of its momentum operator, i.e.. 



^TD 



It can be shown that the EC force [Eq. (pS)) ] follows 
from the TD force expression. Next, we derive the HE 
and IBSC forces that are approximations to the EC force. 



2. The IBS corrected force 

By assuming that the wavefunctions are the eigenstates 
of the Hamiltonian operator, i.e., in matrix form 



(27) 



one can derive an approximation to the EC force which 
we call the incomplete basis set corrected force. Using 
the ground state assumption [Eq. (P7)) ] and the relation 
between the matrices Dq and S, 



dS 
dR^ 



D„ + D! 



(28) 
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one can derive the IBSC force from the EC force [Eq. 

131 



:iIBSC 



dS 



(29) 



Using the above expression, the time derivative of the 
total Hamiltonian reads as 



position-independent basis sets the HE force works 
weU, which is the case, e.g., in the Troulher-Martins 
pseudopotential-based Octopus packagei^. Pseudopoten- 
tials on a position-independent basis also simplify the 
electronic part of the Ehrenfest MD equations [Eq. (|24p ] 
as the P term is zero. In addition to the Octopus pack- 
age, the HE approach is used in the SIESTA package^^ 
which is based on the LCAO basis set. 



IBSC 



ds_ 



-(HS-iD„ + c.c.))c„. 



(30) 

For a linear system, this should be essentially zero be- 
cause the vectors c„ are the eigenstates of the Hamil- 
tonian matrix H. Eor a nonlinear system, however, the 
time derivative is generally non-zero because the differ- 
ence between the gradient of S and the HS"^Da + c.c. 
terms might cause fluctuations in the total energy. Espe- 
cially in the case of nonadiabatic processes, the ground 
state assumption [Eq. (P7| ] might yield inadequate re- 
sults in terms of the total energy conservation. The IBSC 
force has been suggested for Ehrenfest MD in a finite ba- 
sis in literature^. 



3. The Hellmann-Feynman force 

The most traditional way of treating the atomic forces 
is the Hellmann-Eeynman theorem, i.e., the force acting 
on atom a is calculated by taking the expectation value 
of the gradient of the electronic Hamiltonian operator 
with respect to the atom a. Moreover, for the sake of 
consistency (i.e. no wavefunction gradients) , we leave out 
the gradient of the nonlinear term H' in the Hamiltonian 
operator. The HE force is an approximation to the IBSC 
force [Eq. ([^ J as it can be derived by neglecting the 
gradients of the basis functions with respect to the atomic 
positions. Therefore, the HE force in our calculations 
reads as 



riHF 



Now the time derivative of the total Hamiltonian H 
(ITU1)1 can be written as follows 



(31) 
Eq. 



a 71. a 

-5I^--<(HS"'D, + c.c.)c„, (32) 



B. Time propagation of the electron-ion system 

In order to carry out actual simulations, a propagation 
algorithm for the quantum-classical system [Eqs. ((24|) 
and (1^5)) ] is required. Eirst, we use the following splitting 
for the propagation of coupled electrons and ions 

UNjt+At, t) = UN{t + ^,t)U,{t + At, t) 

xUN{t + At,t+^) + OiAt^), (34) 

where the propagator for the nuclei, Uj^, is the stan- 
dard velocity Verlel3, while the electronic states (C/g) 
are propagated using the so-called semi-implicit Crank 
Nicholson (SICN) methodiii^. The idea of the method 
is to first approximate the Hamiltonian matrix H to be 
constant during the time step and solve the following 
linear equation to obtain the predicted future electronic 
states cP""^ 



S + z^(H(0 



cr'(i + Ai) 

c„(t)-)-0(Ai2). (35) 



Then, the predicted future Hamiltonian matrix, based 
on cP''°'^, is used for calculating the Hamiltonian matrix 
in the middle of the time step, H(< -|- At/2) = i(H(t) -f 
TlP'"''^{t+At))+0{At^). Now, using the notation Hi/2 = 
H(i + At/2), the final, propagated states c„(i + At) can 
be obtained from 



S + z^(Hi/2 + P) 



c„(t + At) 



S-^^(Hi/2- 



P) 



c„(t)-HO(At3). (36) 



III. EHRENFEST MD IN A FINITE ELEMENT BASIS 



where the matrix H°'" reads as 



(33) 



Clearly, Eq. shows that the HE force will 

not conserve the total energy because the terms do 
not cancel out each others in general. However, for 



A. The model system 

In order to investigate different forces and their effect 
on the total energy conservation, we study a two-atom 
system using the finite element method (EEM). The ba- 
sis functions depend explicitly on the atomic positions. 
We restrict our calculations to one dimension in order to 
minimize the errors arising from discretization. In order 
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FIG. 1. Non-uniform grid used for the ID model system. The 
grid is shown for two different interatomic distances d = 0.7 
A and d — 1.5 A. For illustratory purposes, the small number 
of FEM basis functions A'' = 50 is used. 



to create a position-dependent basis, we define a map- 
ping / from the uniform grid to a non-uniform grid as 
follows 



fix) = X 



(37) 



where L is the half-width of the grid and k > 1 is a 
parameter. This mapping creates a grid that is more 
dense near the origin and more sparse near the bound- 
aries. Then, we define another mapping g 

z = gifix)) ^f{x) - [{fix) - i?i)e-''(/(-)-«^)' 

+ (/(.t) - i?2)e-''(^'(")-^^)' 

x(l -e-''(^(^)-«i)')], (38) 

where i?i and R2 are the atomic positions, and rjjV > 
are parameters. With the above equation, the points on 
the uniform grid are mapped onto a non-uniform grid 
that is more dense around the atoms. The piece wise 
linear FEM basis functions Xk are then defined on the 
non-uniform grid as follows 



Xk{z) 



Zk + 1 — 



Zk+1- 





if z e [zfe_i,Zfc], 
if z e [zfc,Zfc+i], 
otherwise. 



(39) 



Figure [T] illustrates the non- uniform grid for two different 
interatomic distances, d = 0.7 A and d = 1.5 A. One 
can clearly see that the grid is more dense close to and 
between the atoms, and in turn the sparsity of the grid 
increases as one approaches its boundaries. 

The two-atom system, in our calculations, is defined by 
the Gross-Pitaevskii (GP)-type (see, for example, Ref. [23 
and references therein) energy functional 



(40) 



where p is the electronic density, Ts is the kinetic energy 
of the non-interacting electrons, T4xt is the external po- 
tential and 7 > is a parameter. The external potential 



in our calculations includes the nucleus-nucleus (nn) and 
electron-nucleus (ne) interactions. In order to prevent 
the potential from diverging at the nuclei, we use the 
following, so called soft Coulomb potentials 



Kxt(x) = K: 



nc, soft 

ai 



, soft 



02 



^(x-i?i)2 +ai 
V(i?i -i?2)2+a2' 



(41) 



where oi, 02, ai, 0:2, /3 > are parameters. By varying 
the GP functional [Eq. (j40])] with respect to the density, 
we obtain the following Hamiltonian operator 

H^f + Kc, soft + Kn, soft + 7P- (42) 

Now, the matrix H' [Eq. ([5])] reads as 

H^j - 7 ixMXj) = 7 51 < AuC„, (43) 

n 

where the tensor A is defined as 

X^XJXkXldz. (44) 



Having defined the model system, the forces can be com- 
puted according to the formalism presented in Sec. Ill AI 
and the time propagation proceeds as presented in Sec. 

EH 



B. Results for the model system 

In order to study the different forces presented in Sec. 
Ill AI we present the implementation of the FEM sys- 
tem described in Sec. IIII AI In the example calcula- 
tions, we use the following, reasonable set of parameters 
ai — 1,02 = 3, Q!i = 0.1,0:2 = 0.01,77 = u = 0.8, /? = 
1.2, K = 1.3, and L = 4.0 A. The two lowest states are oc- 
cupied. Moreover, the atomic masses are Mi = 2000 a.u. 
and M2 = 5000 a.u. First, we study a linear system by 
setting 7 = 0, and carry out calculations where the sys- 
tem is initially at the interatomic distance of 1.03 A (the 
equilibrium value being 0.69 A), after which it evolves 
freely for a few femtoseconds. Then, we introduce non- 
linearity to the model system by performing calculations 
with 7 — 0.02 and 7 — 0.2. The total energy conservation 
with the HF, IBSC and EC forces is investigated. Finally, 
we study the validity of the ground state assumption in 
the IBSC force by assigning the kinetic energy of 200 
eV to the atoms which are initially at their equilibrium 
distance. 



1. Linear system 

First, it is instructive to demonstrate the effect of the 
P term on the conservation of the norm of the electronic 
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FIG. 2. Decay of the error in the total energy of the linear ID 
model system as a function of the simulation time step. The 
results obtained with (a) iV = 300 and (b) iV = 50 FEM basis 
functions as well as with the HF, IBSC and EC forces [Eqs. 

, (I29|l and (I25|l , respectively] are shown. The dashed lines 
are quadratic fits, while the blue solid line is just a guide to 
the eye. The results obtained with P = do not fit within 
the scales of the figure. 



states and subsequently total energy conservation. One 
might assume that neglecting this term would merely 
cause small fluctuations in the total energy. In order to 
investigate the validity of this assumption, we simulate 
the two-atom system for T = 3 fs using N = 300 FEM 
basis functions and different time steps. Figure [5] shows 
the maximum fluctuation in the total energy within the 
simulation time, 

AStot = max{|^tot(i) - £^tot(0)| : t G [0,T]}, (45) 

as a function of the simulation time step when the P 
term is not neglected. Clearly, the difference between 
the IBSC and EC forces is negligible as both of them 
result in a quadratic decay of the error in the total energy. 
This result confirms the assumption presented in Sec. 
Ill A 21 that for a linear system the IBSC force should be 
good enough. Moreover, in Fig. [UJa), the total energy 
error with the HF force is roughly as small as that with 
the EC and IBSC forces. This is not surprising as the 
basis is quite large {N — 300). In the case of = 50, 
however, the size of the basis clearly affects the total 
energy conservation as the error is over 100 meV (Fig. 
[2Kb)). Interestingly, in Fig. HKb), the change in the total 
energy for the HF force seems to decrease as the time step 
is increased. However, this is simply due to cancellation 
of errors in this particular case, and for larger time steps 
than shown in Fig. [Ub) the energy change begins to 
increase again. Exactly this behavior is seen in Fig. ^a) , 
but on a different scale. 

By setting P = 0, the behaviour of the total energy 
changes radically as the norm of the electronic states is 
no longer conserved. We find that the maximum fluctu- 
ation in the total energy is 28.3 eV which is 580 % of the 
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FIG. 3. Error in the total energy of the linear ID model 
system with respect to the number of FEM basis functions. 
The results obtained with the HF, IBSC and EC forces [Eqs. 
(|3ip . (|29|) and (|25l) . respectively] are shown. The time step 
used for the calculations is At = 0.05 as. 



maximum kinetic energy of the molecule. Furthermore, 
the norm of the electronic states differs from the value 
of 2 electrons by as much as 0.128 electrons. Thus, the 
P term has a very significant effect on the total energy 
conservation and must not be neglected. 

Finally, we study the validity of the HF theorem by 
comparing the total energy error obtained with the HF 
force to that obtained with the IBSC and EC forces as 
a function of the size of the basis set. These results are 
shown in Fig. [S] We observe that the HF force works suc- 
cificiently well for TV > 100 basis functions. On the other 
hand. Fig. [3] shows that even in a practically adiabatic 
small basis might cause significant fluctuations in 
the total energy. Thus, one has to be careful with the 
number of basis functions in the calculations when us- 
ing the HF force in a position-dependent basis set. In 
order to eliminate the possibility of total energy fluctu- 
ations due to the insufficient number of basis functions, 
one must use either the IBSC or the EC force. 



2. Nonlinear system 

Next, having studied the behaviour of the ID model 
system in the linear case, we turn our attention to the 
nonlinear case by carrying out calculations with a non- 
zero value of the nonlinearity parameter 7 [Eq. (|42[) ]. 

First, similarly to the linear case, we perform calcula- 
tions where the initial atomic separation is d = 1.03 A 
and investigate the total energy conservation. The total 
simulation time is T = 3 fs, and N = 200 FEM basis 
functions are used in all the calculations. Figure |4] shows 
the decay of the error in the total energy with respect 
to the simulation time step. As in the linear case, both 
the IBSC and the EC forces result in a quadratically de- 
caying total energy error with respect to the time step. 
Comparing the results shown in Figs. SJa) andSKb), the 
errors with 7 = 0.02 and 7 = 0.2 appear to be almost 
identical. Furthermore, the HF force works well with 
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FIG. 4. Error in the total energy of the nonlinear ID model 
system as a function of the simulation time step. The results 
obtained with the HF, IBSC and EC forces [Eqs. (gU, ^ 
and (|25[l . respectively] and two different degrees of nonlinear- 
ity 7 = 0.02 (a) and 7 — 0.2 (b) are compared. The dashed 
lines are quadratic fits, while the blue solid line is just a guide 
to the eye. 



N = 200 basis functions as the total energy error due to 
the finite basis is only about 7 meV. 

In order to observe a significant difference between the 
IBSC and EC forces, we have to increase the kinetic en- 
ergy of the system. For this reason, we carry out cal- 
culations where the atoms are initially in equilibrium. 
We then assign the fairly high initial kinetic energy of 
Ek — 200 eV to the atoms, the value of the nonlinearity 
parameter being 7 — 0.2. The total energy conservation 
is investigated within the simulation time T = 1 f s using 
the time step At = 0.02 as. Total energy curves for the 
IBSC and EC forces are presented in Fig. [5] According 
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FIG. 5. Total energy conservation of the nonlinear ID model 
system as a function of the simulation time. The results ob- 
tained with the IBSC and EC forces [Eqs. ^ and ([25|l . 
respectively] are compared. The initial kinetic energy of the 
system is Ek = 200 eV. 



in contrast, causes clear fluctuations in the total energy 
- the error is roughly two orders of magnitudes larger 
than that for the EC force. Moreover, instead of fluctu- 
ating around a constant value, the total energy appears 
to drift downwards as a function of time. Thus, the EC 
force appears to be a much better choice for simulations 
with energetic ions. 



IV. EHRENFEST MD WITHIN THE PAW FORMALISM 

A. Theoretical framework 

Ehrenfest MD within the PAW formalism is similar to 
the finite basis formalism presented in Sec. HIl The first 
notable difference is that within the PAW method, we 
actually have a basis defined by two different functions, 
the projectors p° and the pseudo partial waves 0°, which 
fulfill 



(46) 



where i is a multi-index consisting of the quantum num- 
bers I, m and n. The second notable difference is 
that the dependency on the atomic positions arises from 
the position-dependent PAW transformation operator T- 
Similarly to the finite basis set case, there appears an ad- 
ditional term due to the moving basis set in the TDDFT- 
equivalent of the single-particle TD Schrodinger equation 
[Eq. (in])]. We start from the all-electron TDKS equa- 
tion. 



by applying the PAW transformation. 



(47) 



(48) 



between the all-electron KS wavefunctions tpn and the KS 
pseudo wavefunctions V'n- Then, the all-electron TDKS 
equation [Eq. (|Tf|)] is operated from the left by the ad- 
joint of the PAW transformation operator, T^. Sub- 
sequently, we arrive at the following PAW-transformed 
TDKS equation 



(49) 



where S is the PAW overlap operator, and H is the PAW 
Hamiltonian operator. The P term, which corresponds 
to the P matrix presented in Sec. Ill Al reads as 



(50) 



to the figure, the EC force conserves the total energy very 
well within the simulation time of 1 fs. The IBSC force. 



It takes into account the time evolution of the PAW 
transformation operator in TDDFT-based quantum- 
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classical MD simulations. Qian et al.— derived the fol- 
lowing expression for this term 

*5]v„.(i + £t)_^(i + 4) 

a ^ 

=-j^Va-D„ (51) 

a 

where ia = X^id'^i') ~ I'i'i)) iPil ^ projection operator 
belonging to atom a, and we have defined the operator 
Da in the spirit of the formalism presented in Sec. HIl 
Moreover, Eq. ((5T|) only holds if the overlap between the 
PAW augmentation spheres is zero. In practice, how- 
ever, Eq. ([?T|) turns out to work well even in the case of 
overlapping augmentation spheres as long the overlap is 
not significant. The operator can be written in the 
following form [Appendix \X\ 

il,i2 " 

d(hf ~ dd>f 

+\pvm^^j-{<^t\^j){pi\]- (52) 

The matrix elements describe the overlap between 

the all-electron and pseudo partial waves 

o^,,.. = (CIC)-(CIC)- (53) 

The next task is to derive a computable expression for 
the force similar to the EC force [Eq. (P5|)]. Using the 
same reasoning as in Sec. |lTl i.e., the conservation of the 
total energy, we replace the IBSC force expression used 
for ground state calculations in the GPAW package22i2i, 

pIBSC = -^^+Y: (^«I^IV^.) , (54) 

where /„ is the occupation number of state n, with the 
more general expression 

F^^ = + E /" i^n\^VS-'H + c.c.\i^n) . (55) 

n 

By defining cjn — S^^Hipn and the vector- valued matrix 
elements 

AL.. = (Clg|;)-(CJ^), (56) 
we arrive at the following equation 

+ {^n\pl) Al,Jpl\~9n)]. (57) 

With the EC force expression for GPAW and the re- 
quired force corrections [Eqs. (1551) and ([57]) ] . the atomic 
forces can be straightforwardly calculated. The matrix 



elements Of , and Af , are calculated on radial grids 
inside the PAW augmentation spheres, whereas the terms 
involving the pseudo wavefunctions are computed on uni- 
form Cartesian grids. 

The coupled electron-ion system is propagated using 
the SICN method [Eqs. and for the electronic 
states ipn and the velocity Verlet algorithm for the nuclei 
in a similar fashion as in the ID calculations in the FEM 
basis. However, unlike in one dimension, the calculation 
of the inverse PAW overlap operator is not trivial. 
We use two different approaches for this. In the first 
approach we assume that the inverse overlap operator 
can be written as 

S-' = 1 + Y.T.\pVCI,ApI\' (58) 

a ii,i2 

where C," are the inverse overlap coefficients. They 
can be obtained from the linear equation resulting from 
the requirement 

S-^S^l, (59) 

in which we also assume that there is no overlap between 
the PAW augmentation spheres. For this reason, this 
approach will likely produce unsatisfactory results if the 
overlap is non-zero. In the second approach, we obtain 
the required S^^Htpn terms by solving the linear equa- 
tions 

Sx = Hi^n, (60) 

iteratively using the conjugate gradient (CG) method. 
This approach is generally more accurate than the app- 
roximative inverse method [Eqs. (f58| and (f59| ] and less 
sensitive to PAW overlap and numerical errors due to 
a large grid spacing. This is illustrated by Table U in 
which the approximative inverse method is very sensitive 
to the grid spacing used in the calculations. However, 
being an iterative approach, the CG method will likely 
consume more computational resources than the approx- 
imative method. 



B. Calculations for small and medium-sized molecules 

In order to test our PAW based Ehrenfest MD method, 
we study the dynamics of small and medium sized 
molecules both in adiabatic and in nonadiabatic cases. 
The adiabatic cases include the vibration of the NaCl 
molecule, and the rotation of the H2C = NH^ molecule 
about its internal axis with a small initial kinetic energy. 
Nonadiabatic effects, then, are studied in the case of H2C 
= NHJ with a high initial kinetic energy, and also in the 
hydrogen bombardment of the C40H16 molecule. We use 
the LDA exchange-correlation functional2^ in all the cal- 
culations. The grid spacing is /i = 0.2 A unless specified 
otherwise. 
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1. Vibration of the NaCI molecule 

First, we study the total energy conservation of a sim- 
ple dimer, the NaCl molecule. The simulations begin 
from equilibrium (doq = 2.36 A) with the initial kinetic 
energy of 2 eV. The period of the NaCl vibration calcu- 
lated with the BOMD is Tbomd = 181.6 fs. With this 
low initial kinetic energy, the vibration is almost adi- 
abatic, and subsequently the period obtained with the 
Ehrenfest MD (At = 8 as, IBSC force) is very close to 
the BOMD resuh, Tef = 181.4 fs. This is not surpris- 
ing as the Ehrenfest MD should reduce to the BOMD 
in the case of adiabatic processes. We investigate the 
difference between the IBSC and EC forces in terms of 
the total energy conservation, using both the approxi- 
mative and the CG inverse method for calculating the 
inverse overlap operator. The decay of the error in the 
total energy as a function of the simulation time step is 
shown in Fig. [6] For time steps below 10 as, there is 



gies, Ek = 1.5 eV and Ek = 10 eV. In order to investigate 
the nonadiabaticity in our simulations, we also carry out 
PAW-based BOMD calculations for both initial kinetic 
energies. Based on the results presented in Ref. [2^, the 
rotation is expected to be adiabatic with Ek = 1.5 eV, 
whereas with Ek ~ 10 eV we expect the Ehrenfest MD 
results to clearly differ from the BOMD results. The 
molecule is initially in the planar equilibrium geometry. 
For both initial kinetic energies, we carry out calculations 
using two different time steps. At — 2 and 5 as. Similarly 
to the calculations for the NaCl molecule, the IBSC and 
the EC forces are applied. In the case of the EC force, 
both the approximative and the CG method are used for 
computing the inverse overlap operator . 

The potential energy surfaces obtained from the 
Ehrenfest MD and the BOMD simulations for both ini- 
tial kinetic energies are presented in Fig. [T] For the 
Ehrenfest PES, the EC force in conjunction with the CG 
method for calculating S^^ is used. The figure illustrates 
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FIG. 6. Decay of the error in the total energy of the NaCl 
molecule as a function of the simulation time step. The re- 
sults are shown for the IBSC [Eq. ^] and the EC [Eq. (f55|l ] 
forces. Approx denotes the approximative method for calcu- 
lating the inverse overlap operator [Eqs. (|58[) and H59p ]. and 
CG corresponds to the conjugate gradient method. The lines 
are just a guide to the eye. 

little difference between the IBSC and EC forces. Actu- 
ally the total energy error with the IBSC force is slightly 
smaller than that with the EC force, thus rendering the 
EC force unnecessary in this adiabatic case. Moreover, 
the two methods used for calculating the inverse overlap 
operator give pratically identical results. 



2. Rotation of the H2C — NHt molecule 



Next, we turn our attention to simulating nonadia- 
batic dynamics. The dynamics of the H2C — NH^ 
molecule has been studied with both the Hartree-Fock 
based Ehrenfest dynamics^^ and the trajectory surface 
hopping method^ . We study the rotation of this molecule 
about its internal axis by carrying out Ehrenfest MD cal- 
culations with two different initial torsional kinetic ener- 
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FIG. 7. Time evolution of the PES of the H2C = NHJ 
molecule. The results obtained with the Ehrenfest MD and 
the BOMD for low and high initial kinetic energies are com- 
pared. In the BOMD simulations, the time step of 0.1 fs is 
used. The Ehrenfest MD simulations are performed using the 
EC force [Eq. (|55p ] in conjunction with the CG method for 
calculating the inverse overlap operator. 

that the dynamics is nearly adiabatic with Ek = 1.5 eV 
as the Ehrenfest and the Born-Oppenheimer PES are al- 
most identical. In contrast, with Ek = 10 eV, the Ehren- 
fest PES starts to deviate rapidly from the BO one, which 
indicates a significant amount of nonadiabaticity in the 
dynamics. 

In the adiabatic case, the total energy is conserved to 
a few meV even with the IBSC force. With Ek — 10 
eV, in contrast, this is no longer the case. The total 
energy curves obtained with the different forces and the 
time steps of 2 and 5 as are shown in Fig. [51 First, 
we observe that the maximum total energy fiuctuation 
with the IBSC force is roughly the same for both time 
steps used in the calculations. Secondly, the EC force 
in conjunction with the conjugate gradient method for 
calculating the inverse overlap operator, yields the best 
total energy conservation for both time steps, Ai?tot = 
45 meV and 9.3 meV for At = 5 and 2 as, respectively. 
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FIG. 8. Conservation of the total energy of the H2C = NHj 
molecule as a function of the simulation time. The results 
obtained with the IBSC and the EC forces [Eqs. and 
(I55|l . respectively] using two different time steps, At — 2 and 5 
as, are compared. The initial kinetic energy of the molecule is 
Ek = 10 eV. Approx and CG correspond to the approximative 
[Eqs. (|58|l and (|59[) ] and the CG method for calculating the 
inverse overlap operator, respectively. 



10-20 % more computational time, the CG inverse should 
be used for simulations involving nonadiabatic effects ins- 
tead of the approximative inverse. 



3. Hydrogen bombardment of the C40H16 molecule 

As the final test for our PAW-based Ehrenfest MD 
method, we study the collision of hydrogen with the 
graphene-like nanoflake C40H16. With high enough im- 
pact energies, the hydrogen projectile will invoke elec- 
tronic excitations in the target, rendering the BOMD 
approach unusable. In Ref. |9|, a similar hydrogen ion 
stopping process was studied in the case of graphene, for 
two representative trajectories, 1) center of hexagon and 
2) the impact parameter of 0.25 A along the C-C bond. 
The two trajectories are illustrated in Fig. |9l The EC 




The latter number is very good considering the amount 
of nonadiabaticity involved in the dynamics, and it could 
be further improved by decreasing the time step. Fur- 
thermore, with At = 2 as, the EC force in conjunction 
with the approximative inverse already improves the to- 
tal energy conservation quite significantly compared to 
the IBSC force. Nevertheless, the error in the total ener- 
gy is at least twice as high as with the CG inverse. 

Next, we study the effect of the grid spacing on the to- 
tal energy conservation by carrying out simulations with 
h = 0.15, 0.2, 0.25 A. Because the IBSC force works well 
with the low initial kinetic energy of Ek = 1.5 eV, we only 
study the more energetic case. The error in the total en- 
ergy as a function of the grid spacing is presented in Table 
|T1 The accuracy of the approximative inverse increases 

TABLE I. Maximum total energy fluctuation of the H2C = 
NH^ molecule as function of the grid spacing h. The initial 
kinetic energy is Ek ~ 10 eV. The results obtained with the 
IBSC and EC forces [Eqs. (IMl) and ((SSJ, respectively] are 
compared. All the energies are in meV. 
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149.19 
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151.06 


228.78 
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rapidly as a function of decreasing grid spacing. With 
h = 0.25 A, the IBSC force actually conserves the total 
energy better than the EC force if the approximative in- 
verse is used. However, despite the rapid convergence, 
it is undesireable that the total energy error has such 
a strong dependence on the grid spacing. In contrast, 
the grid spacing has very little effect on the results ob- 
tained with the CG inverse. For this reason, even though 
the calculations might require, depending on the system. 



25 ^ 



FIG. 9. C40H16 molecule and the two representative trajec- 
tories used in the simulations. 

force in conjunction with the CG method for computing 
the inverse overlap operator is used in all the calcula- 
tions. The time step varies between 1 and 5 as such that 
At = 1 as is used in the simulation with the highest ini- 
tial projectile energy, Ek = 10 keV and At = 5 as in that 
with the lowest initial energy, Ek — 100 eV. 

We study the accomodation of the energy transferred 
into the individual degrees of freedom. Because of the 
overlap between the augmentation spheres, the PAW 
method underestimates the atomic forces when the in- 
teratomic distance is small. Consequently, we use pair- 
potential corrections derived from results obtained with 
FHl-aims2I when the distance between the projectile and 
the nearest carbon atom is smaller than 0.5 A. Figure 
[TUf a) shows the total transferred energy and the C re- 
coil energy for the bond trajectory. Electronic excita- 
tions significantly influence the results beyond the im- 
pact energy of 400 eV as the total transferred energy 
starts to increase. This observation is in agreement 
with the Troullier-Martins pseudopotential calculations 
for graphene in Ref. H. However, despite the qualitative 
agreement, the quantitative results differ slightly, which 
can be attributed to the heavy overlap between the aug- 
mentation spheres of the hydrogen projectile and the tar- 
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get carbon atom. Nevertheless, the agreement is quite 
good considering that the Eqs. ([5T|) and ([5^ used for 
calculating the P operator are in principle only correct 
when the augmentation spheres do not overlap. 




100 1000 100 1000 10000 



H atom initial energy (eV) H atom initial energy (eV) 

FIG. 10. Deposition of the energy into the individual degrees 
of freedom in the collision of the H atom with the C40H16 
molecule, (a) Energy transferred to the target C atom, and 
the total transferred energy as a function of the H impact en- 
ergy. The results are shown for the bond trajectory. Pair-p 
denotes the pair-potential corrections, (b) Energy transferred 
into the electronic degrees of freedom as a function of H im- 
pact energy for both trajectories. The hues are just a guide 
to the eye. 



V. CONCLUSIONS 

We have described the implementation of the Ehren- 
fest molecular dynamics within the time-dependent den- 
sity functional theory and the projector augmented-wave 
method. Moreover, we have studied different Ehrenfest 
MD forces using a position-dependent finite element basis 
in one dimension as well as in 3D with the GPAW pack- 
age. In the ID calculations, the forces were compared by 
carrying out Ehrenfest MD simulations for a two-atom 
system. In the 3D calculations, the dynamics of small 
and medium sized molecules was studied both in adia- 
batic and in nonadiabatic cases. The incomplete basis 
set corrected and Hellmann-Feynman forces were found 
to work succifiently well in adiabatic or nearly adiabatic 
cases, whereas in clearly nonadiabatic cases unphysical 
fluctuations in the total energy were observed. The to- 
tal energy conserving force was found to work well also in 
the nonadiabatic regime. Finally, the PAW-based Ehren- 
fest MD results were compared to TrouUier-Martins pseu- 
dopotential calculations. From these results, we conclude 
that our method seems applicable to simulating the nona- 
diabatic dynamics of medium sized molecules and beyond 
as long as the PAW augmentation spheres do not overlap 
significantly. 



In Fig. [rUrb). the energy transferred into the electronic 
degrees of freedom is presented. The PAW and TM re- 
sults are in good agreement for the center of hexagon tra- 
jectory, i.e., when there is no overlap between the PAW 
augmentation spheres. Thus, our method describes elec- 
tronic excitations in a similar fashion as the TM based 
method. Consequently, it seems that our method works 
correctly also in nonadiabatic cases as long as the over- 
lap between the PAW augmentation spheres is not sig- 
nificant. Even in the case of overlapping augmentation 
spheres, our method can predict qualitative trends with 
the help of pair-potential corrections. 

Finally, we summarize the results regarding the total 
energy conservation in the calculations. First, in all the 
simulations for the center of hexagon trajectory, the total 
energy is conserved to better than 2.7 meV, which is an 
excellent result considering that the amount of energy 
deposited into the electronic degrees of freedom is of the 
order of tens of eV. In the case of the bond trajectory, 
unfortunately, such good numbers cannot be obtained 
- the total energy is conserved to better than 230 meV. 
This is probably due to the breakdown of the zero overlap 
approximation in deriving the P term as it is essential 
that this term is correct in order to conserve the total 
energy. Nevertheless, significant PAW overlap is quite 
unusual in Ehrenfest MD applications. Thus, in most 
cases our method can be expected to conserve the total 
energy very well. 
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Appendix A: Symmetric form of the operator Da 

In order to carry out Ehrenfest MD simulations in 
practice, it is useful to write the P term [Eq. ([5T|) ]. in 
a symmetric form reminiscent of the other observables 
within the PAW method. This can be achieved by de- 
riving a symmetric form for its constituent operator 
[Eq. ([5^ ]. We start by expanding this operator in terms 
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of the PAW projectors and partial waves 



rif r) 



E IK) ((CI -(CI) 



symmetric form of the operator, 



dRa 



E(ic)-ic))(c 



(Al) 



Rearranging the terms and adding and substracting a 
suitable term gives 



= (1 - E m (CD^R^EdC) - ic)) iPi 



■Eic)(ci7^E(ic)-ic))(p-j 



-EiC)(Ci7m-E(iC)-iC))(?5?. 



dRa 

d_ 
dRa 



EiPn)(Cia|;E(iC)-iC))(p-, 



(A2) 



Using the orthonormality of the projectors and pseudo 
partial waves, we obtain the following expression 



D. = -ElK)(ClaR;ElC)(C 



ElC)(ClaR;ElC)(C 



EiC)(Ci^Eic)(c 



+ E w (C I E(ic) - ic)) (p". I • (A3) 



The first term in Eq. (|A3|) is zero due to orthonormality. 
Combining the remaining terms then yields the desired 
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